A new approach to the complex-action problem 
and its application to a nonperturbative study of superstring theory 
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Monte Carlo simulations of a system whose action has an imaginary part are considered to be 
extremely difficult. We propose a new approach to this 'complex- act ion problem', which utilizes 
a factorization property of distribution functions. The basic idea is quite general, and it removes 
the so-called overlap problem completely. Here we apply the method to a nonperturbative study of 
superstring theory using its matrix formulation. In this particular example, the distribution function 
turns out to be positive definite, which allows us to reduce the problem even further. Our numerical 
results suggest an intuitive explanation for the dynamical generation of 4d space-time. 
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I. INTRODUCTION 

It occurs in many interesting systems ranging from 
condensed matter physics to high-energy physics that 
their action has an imaginary part. Some examples for 
instance in high-energy physics are the finite density 
QCD, Chern-Simons theories, systems with topological 
terms (like the 0-term in QCD), and systems with chi- 
ral fermions. While this is not a conceptual problem, it 
poses a technical problem when one attempts to study 
these systems by Monte Carlo simulations, which would 
otherwise provide a powerful tool to understand their 
properties from first principles (see Refs. 0,ai for re- 
cent works). 

In this Letter we propose a new approach to this 
'complex-action problem'. Suppose we want to obtain 
an expectation value of some observable. Then, as a 
more fundamental object, we consider the distribution 
function associated with that observable. In general the 
distribution function has a factorization property, which 
relates it to the distribution function associated with the 
same observable but calculated omitting the imaginary 
part of the action. The effect of the imaginary part is 
represented by a correction factor which can be obtained 
by a constrained Monte Carlo simulation. One of the 
virtues of this method is that it removes the so-called 
overlap problem completely. This problem comes from 
the fact that the two distribution functions — one for 
the full model and the other for the model omitting the 
imaginary part — have little overlap in general. The 
method avoids this problem by 'forcing' the simulation 
to sample the important region for the full model. 

The determination of the correction factor becomes in- 



creasingly difficult as the system size increases. In this 
sense, our approach does not solve the complex action 
problem completely. This should be contrasted to the 
meron-cluster algorithm p], with which one can study a 
special class of complex-action systems by computer ef- 
forts increasing at most by some power of the system size. 
The factorization method eliminates the overlap prob- 
lem, which composes some portion of the complex action 
problem, but not the whole. However, the resolution of 
the overlap problem is in fact a substantial progress. For 
instance, Refs. |2j developed a new method to weaken 
the same problem in finite density QCD, and the criti- 
cal point was successfully identified. Therefore we expect 
that the complete resolution of the overlap problem al- 
lows us to address various interesting questions related 
to complex-action systems with the present computer re- 
sources. Since our method is based on the general prop- 
erty of distribution functions, it can be applied to any 
complex-action systems. 

In this article we are concerned with a nonperturbative 
study of superstring theory using its matrix formulation 
|4j. Eventually we would like to examine the possibility 
that our 4-dimensional space time appears dynamically 
in 10-dimensional string theory @, 0, Q, S H Il3 • Monte 
Carlo simulation of the matrix model suffers from the 
complex action problem, and there are evidences that 
the imaginary part of the action plays a crucial role in 
the dynamical reduction of the space-time dimensionality 
0- We will discuss how we can study such an issue by 
Monte Carlo simulation using the new approach. 



II. 



THE SUPERSTRING MATRIX MODEL 
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As a nonperturbative definition of type IIB superstring 
theory in 10 dimensions, Ishibashi, Kawai, Kitazawa and 
Tsuchiya Q proposed a matrix model, which can be 
formally obtained by the zero- volume limit of D = 10, 
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J\f = 1, pure super Yang-Mills theory. The partition 
function of the type IIB matrix model (and its obvious 
generalizations to D = 4, 6) can be written as 



Z = dA 



Zt[A] , 



(1) 



where (fi 



1, • • • , D) are D bosonic N x N traceless 
hermitian matrices, and Sb = — ip'Tr (L4 M , A u ] 2 ) is the 
bosonic part of the action. The factor Zf [A] represents 
the quantity obtained by integration over the fermionic 
matrices, and its explicit form is given for example in 
Refs. 0, El- The convergence of the integral for 
arbitrary N > 2 was first conjectured 0] an d proved 
recently [TsT ] . The only parameter g in the model can be 
absorbed by rescaling A^ i— > y/gA^ , which means that g is 
merely a scale parameter rather than a coupling constant. 
Therefore, one can determine the q d ependence of any 
quantities on dimensional grounds [la] • Throughout this 
paper, we make our statements in such a way that they 
do not depend on the choice of g. 

In this model space-time is represented by A^, and 
hence treated dynamically p| • It is Euclidean as a result 
of the Wick rotation, which is always necessary in path 
integral formalisms. Its dimensionality is dynamically 
determined and can be probed by the moment of inertia 
tensor defined by 



V = lTr(V4„) 



(2) 



Since T^ v is a D x D real symmetric matrix, it has D 
real eigenvalues corresponding to the principal moments 
of inertia, which we denote as Ai (i = 1, • • • , D) with the 
ordering 



Ai > A 2 > ■ • ■ > Xd > 



(3) 



Let us define the VEV (O) with respect to the partition 
function If we find that (Ai) with i — l,---,d is 

much larger than the others, we may conclude that the 
dimensionality of the dynamical space-time is d . 



III. THE COMPLEX ACTION PROBLEM 

The fermion integral Zf[A] in the partition function 
Q is complex in general for D = 10, N > 4 and for 
D = 6, N > 3 Q. Let us restrict ourselves to these cases 
in what follows. Parameterizing the fermion integral as 
Zf{A] = exp(rR + «r), the partition function may be 
written as 



dA e 



-So e ff 



(4) 



where So = Sb — Tr is real. According to the standard 
reweighting method, one evaluates the VEV (Ai) as 



(K) = 



(5) 



where the symbol ( • )o denotes a VEV with respect to 
the partition function 



Zo 



dA 



-So 



(6) 



The VEV ( • )o can be evaluated by standard Monte Carlo 
simulations. However, (e* r ) is nothing but the ratio of 
the two partition functions Zq/Z, and therefore it be- 
haves as AF at large N, where AF > is the differ- 
ence of the free energy density of the corresponding two 
systems. This enormous cancellation (note that |e* r | = 1 
for each configuration) is caused by the fluctuation of 
the phase T, which grows linearly with the number of 
fermionic degrees of freedom, which is of 0(N 2 ). As a 
result the number of configurations required to obtain the 
VEV (e* r ) with sufficient accuracy grows as e constJV . 
The same is true for the numerator ^Ai e ir ) in JSJ. This 
is the notorious 'complex action problem' (or rather the 
'sign problem', as we see below), which occurs also in 
many other interesting systems. 

In fact we may simplify the expression (JSJ by using a 
symmetry. We note that under parity transformation : 



A? 



-Ax 



Af = A, for 2 < i < D , 



(7) 



the fermion integral Zf[A] becomes complex conjugate 
while the bosonic action Sb is invariant. Since the 
observable Ai is also invariant, we can rewrite (JSJl as 



(A, 



(Ai cosT) 
(cosr) 



(8) 



Note, however, that the problem still remains, since cosT 
flips its sign violently as a function of A^. 



IV. THE NEW METHOD 
A. The factorization property of distributions 

The model (0 omitting the phase T was studied up to 
N = 768 and N = 512 for D = 6 and D = 10 respec- 
tively using the low-energy effective theory |jg. There 
it was found that (A^o/^iV 1 / 2 ) approaches a universal 
constant independent of i as iV increases. This means 
that the dynamical space-time becomes isotropic in D- 
dimensions at iV = oo, and hence the absence of SSB of 
SO (I?) symmetry, if one omits the phase T. 

We normalize the principal moments of inertia Ai as 



A,; — 



A, 



(A 



i/0 



(9) 



Then the deviation of (Ai) from 1 represents the effects of 
the phase. The relevant question is whether the deviation 
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depends on i at large N. In order to obtain the expec- 
tation value (Aj), we consider the distribution associated 
with the observable Aj : 



Pi[x) = f (S(x ~ Aj)) 



(10) 



As an important property of the distribution pi(x), it 
factorizes as 



Pi(x) = q Pp (x) Wi(x) , 
where C is a normalization constant given by 

C= <e ir ) = (cosr) . 
The real positive function pf^ (x) is defined by 

pf(x)^(^-Aj)) , 



(11) 



(12) 



(13) 



which is nothing but the distribution of Aj in the model 
© without r. The function p^°\x) is peaked at x = 1 
due to the chosen normalization JUJ- The function Wi(x) 
in (|ll|l can be regarded as the correction factor repre- 
senting the effect of T, and it is given explicitly as 



Wi{x) d = (e ir )j, x = (cosT)i tX , 



(14) 



where the symbol ( • )i^ x denotes a VEV with respect to 
a yet another partition function 



Z ifX = I dAe~ So S(x - A,) 



(15) 



Given all these definitions, it is straightforward to prove 
the relation (|llfl . 



B. Monte Carlo evaluation of pf*\x) and wt(x) 

In order to obtain the function Wi (x) , we have to sim- 
ulate (|15fl . In practice we simulate instead the system 



Z it v =[dA e" s ° e~ v ^ , 



(16) 



where V(Xi) is some potential introduced only for the 
i-th principal moment of inertia. The explicit form of 
the potential we used in the study is V(z) = \l{ z — 
£) 2 , where 7 and £ are real parameters. The results are 
insensitive to 7 as far as it is sufficiently large and we took 
7 = 1000.0. Let us denote the VEV associated with the 
partition function (|16fl as (O)iy. Then the expectation 
value (cosT^y provides the value of Wi (x) at x = (\i)iy. 

The function pf 1 ' (x) can be obtained from the same 
simulation l|16fl . Note that the distribution function for 
Xi in the system <|16[) is given by 

Pi,v(x) d = (Six-Xi))^ cx pf'(i)e- v (W' a ) . (17) 



The position of the peak x p is given by the solution to 
= JLlnp iy (x) = f^ix) - (Ai)oV"((A<)oa:) , (18) 



where we have defined 

f(0) 



}r^^-\npf\x) 



(19) 



This implies that (Xi)aV {{Xi)oX v ) gives the value of 

fi(x) at x = x p . Since we take 7 sufficiently large, 
the distribution piyix) has a sharp peak, and we can 
safely replace the position of the peak x p by the expecta- 
tion value (Xi)i.v- Once we obtain (x), we can obtain 
pf\x) by 



(0)/ ^ 
PI [x) = exp 



/ dz ff" 1 (z) + const. 
Jo 



(20) 



where the integration constant can be determined by the 
normalization of p[°\x). 



C. Resolution of the overlap problem 



by 



From pf*\x) and Wi(x), we may obtain the VEV (Aj 



(Aj 



J Q dxp\ '{x)wi(x) 

Actually this simply amounts to using the reweighting 
formula © but calculating the VEVs on the r.h.s. by 



(Aj cos r) 



(cosr)o 



dxxpf\x)wi{x) (22) 



dx p\ ^ (x) w % (x) 



(23) 



This reveals one of the virtues of our approach as com- 
pared with the standard reweighting method using the 
formula (JSJ) directly. Suppose we are to obtain the l.h.s. 
of and (|2*3*|l by simulating the system ©. Then 
for most of the time, Aj takes the value at the peak of 
Pi (x). However, in order to obtain the VEVs accurately 
we have to sample configurations whose Aj takes a value 
where \p[°\x)wi(x)\ becomes large. In general the over- 
lap of the two functions becomes exponentially small with 
the system size, and this makes the important sampling 
ineffective. Therefore, this 'overlap problem' composes 
some portion of the complex-action problem. The new 
approach eliminates this problem by 'forcing' the simu- 
lation to sample the important region. 



D. Further improvement in the case Wi(x) > 

So far, we have been discussing the general properties 
of the new method. In the case at hand, we can actually 
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FIG. 1: The function In Wi(x) is plotted for TV =12,16,18,20. For x < 1 we also plot data for N = 4, 8 to clarify the 
convergence. We extract the scaling function $4(2) by fitting the data to some analytic function, which is represented by the 
solid line. The dashed line represents "l>5(x), which is obtained similarly from the scaling behavior of In W5 (x). 



further reduce the problem by using the fact that the 
correction factors Wi (x) are actually positive definite, and 
so is the full distribution function Pi(x). (Note that this 
is not guaranteed in general.) This allows us to obtain 
the VEV (Xi) by minimizing the 'free energy density' 
Fi(x) = —-j^i log pi (x), instead of using For that 

we simply need to solve F[(x) — 0, which is equivalent to 



N* h [ ' dx 



(24) 



The function in the bracket [ • ] is expected to approach 
a well-defined function as N increases : 



1 

N2 



lnwi(x) — > $i(x) 



(25) 



Let us note that Wi(x) is nothing but the expectation 
value of e lT in the system (11511 . According to the ar- 
gument below ©, Wi(x) for fixed x decreases as e~ aN 
at large N. The constant a may depend on x, hence 
the assertion. Indeed our numerical results in Fig. ^ (al- 
though the achieved values of N are not very large) seem 
to support this argument. Once we extract the scaling 
function &i(x), we may use it instead of ■^j\nwi(x) in 
(|2~4"|| for larger N. Thus we are able to obtain the VEV 
(Xi) for much larger V than those allowing the direct 
Monte Carlo evaluation of the correction factor Wi(x). 

The positive definiteness of Wi (x) is crucial for such an 
extrapolation technique to work. If we were to calculate 
the VEV (Xi) by (|21|l . we would need to calculate the 



correction factor for larger N by Wi(x) = e N ®*( x \ where 
the multiplication by N 2 and the exponentiation would 
magnify the errors in <frj(x) considerably. This does not 
occur when we obtain the VEV (Xi) by solving (|24|l . 



V. RESULTS 

Monte Carlo simulation of Ijlrjfl can be performed by 
using the algorithm developed for the model © in Ref. 
[Tl|. The required computational effort is 0(iV 6 ). In 
this work, we use instead the low-energy effective thcor 
proposed in Ref. [l| and further developed in Ref. 
The required computational effort becomes 0(/V 3 ). For 
the definition of the low-energy effective theory as well 
as all the technical details including parameters used in 
the simulations, we follow Ref. The validity of the 
low-energy effective theory in studying the extent of the 
dynamical space time is discussed in Ref. 0. We also 
note that the complex-action problem survives in passing 
from the full theory to the low-energy effective theory, 
and hence we expect that the effects of the phase on the 
reduction of space-time dimensionality should be visible 
also in the low-energy effective theory, if it is there at all. 
Here we study the D = 6 case (instead of D — 10, which 
corresponds to the type IIB matrix model) to decrease 
the computational efforts further. 

In Fig. we plot -^2 lnw 4 (x). The correction factor 
Wi(x) has a minimum at x ~ 1 and it becomes larger for 
both x < 1 and x > 1. This can be understood as follows. 
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FIG. 2: The function ^ff\x) is plotted for N = 64, 128. 
scaling function $4(2:) extracted in Fig. 



The solid line represents —$4 (a;), which we calculate from the 



Let us recall again that Wi (x) is the expectation value of 
e ,r in the system (|15H , where A^ is constrained to a given 
value of x. At x — 1, the system Ijl5|) is almost equivalent 
to the system ©, because Ai would be close to 1 even 
without the constraint. (From this, it also follows that 
Wi(l) takes almost the same value for all i.) Therefore, 
the dominant configurations of the system Ijl5(l at x = 
1 is isotropic at large N On the other hand, the 
dominant configurations of the system i|15|) at small x 
are (i — l)-dimensional, since the constraint forces A^ to 
be small, and due to the ordering © , all the Xj with j > 
i become small. Similarly the dominant configurations 
of the system (|15f) at large x are almost i-dimensional, 
since the constraint forces Ai to be large, and due to the 
ordering (J3J, all the Xj with j <i become large. Now let 
us recall that the phase T vanishes when the configuration 
A has the dimensionality d < <i cr , where d CI = 4, 6 for 
D — 6, 10, respectively 0- As a consequence, Wi(x) gets 
larger in both x < 1 and x > 1 regimes. 

As mentioned already, Fig. ^ supports the scaling be- 
havior (|25[) with increasing N. The scaling function 
$4(2:) can be extracted by fitting the data to some ana- 
lytic function. We find that $4(2;) approaches linearly 
as x — > 0, and it approaches some negative constant expo- 
nentially as x — > 00. We observe a similar scaling behav- 
ior for In ws(x). The corresponding scaling function 
$5(2;) is plotted in Fig. for comparison. 

Fig. [5] represents a graphical solution of 124|) for i = 
4. The open and closed circles describe the function 
777/4°^ (x) for N — 64,128 respectively. It is positive 
at x < 1 and turns negative at x > 1, which reflects 



the fact that p < f\x) is peaked at x = 1. The solid 
line represents — $ 4 (x). The intersections of the two 
curves provide the solutions to (|24|l . At N — 128, we 
find that the distribution p4(x) has two peaks; one at 
x — x s < 1 and the other at x = xi > 1. The ratio 
of the peak height R = pi{x s ) / p^{xi) can be written as 
R = cxpjiV 2 (A s — Ai)}, where A s and Ai are the area 
of the regions surrounded by the two curves. We obtain 
A s ~ 5.0 x 10~ 4 and Ai ~ 4.5 x 10~ 3 , from which we 
conclude that the peak at x > 1 is dominant. In Fig. 
we show the results of a similar analysis for p$(x). We 
find that the distribution p$(x) at N = 128 also has two 
peaks; one at x < 1 and the other at a; > 1. However, 
here we obtain A s - 2.0 x 10~ 3 and Ai ~ 3.8 x 10~ 3 , 
which are comparable. 



VI. SUMMARY AND DISCUSSIONS 

We have proposed a new method to study complex- 
action systems by Monte Carlo simulations. In partic- 
ular we discussed how we can use the method to inves- 
tigate the possibility that four-dimensional space time 
is dynamically generated in the type IIB matrix model. 
The space-time dimensionality is probed by the eigenval- 
ues Ai of the moment of inertia tensor and we study the 
distribution of each eigenvalue. The distribution pf*\x) 
obtained without the phase T has a single peak, which 
is located at x = 1. The effect of the phase T on the 
distribution function is represented by the multiplication 
of the correction factor Wi(x) as stated in (|llfl . 
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FIG. 3: The function ^/ 6 (0) (a;) is plotted for N = 64, 128. 
scaling function "3>s(x) shown in Fig. 



The solid line represents —& s (x), which we calculate from the 



Our results for the 4-th and 5-th eigenvalues (i = 4, 5) 
in the D = 6 case show that the correction factor Wi(x) 

strongly suppresses the peak of pf\x) at x = 1 and 
favours both smaller x and larger x. As a result, we ob- 
serve that the distribution pt{x) including the effects of 
the phase, in fact, has a double peak structure. More- 
over, the two peaks tend to move away from x = 1 as N 
is increased. It is important to determine which of the 
two peaks becomes dominant in the large N limit. At 
N = 128, we observe that the peak at x > 1 is dominant 
for both pi{x) and ps(x). We note, however, that it is 
much more dominant for pa{x) than for p§(x). 

From Figs. |21 and [31 we observe that the function 
1 fi°\x) changes drastically as we go from N = 64 



N 2 

to N 



128. In fact we find that jjf> ' (x) scales (notice 
the normalization factor j^), as shown in Fig.Q]for i = 4. 
The scaling region extends from x ~ 1, where jjfi {x) 

crosses zero, namely the place where p,- "* (x) has a peak. 
A similar scaling behavior is observed for i = 5. This 
scaling behavior is understandable if we recall that the 
long-distance property of the system is controlled by a 
branched-polymer like system |5j], which is essentially a 
system with N degrees of freedom. If we naively extrap- 
olate this scaling behavior of jrf^(x) to larger N, the 
l.h.s. of (|24|l becomes negligible. It follows that the peak 
at x < 1 eventually dominates for both i = 4, 5, consid- 
ering the asymptotic behaviors of $i(x) as x — > and 
x — > oo. This means that the space-time dimensional- 
ity becomes d < 3. However, it is well-known that the 
Hausdorff dimension of a branched polymer is c?h =4, 



which implies that such a system is not easy to collapse 
into a configuration with dimensions < 3. The conse- 
quence would be that pf^ (x) is much more suppressed in 

the small x regime than p^ (x) at large N. We consider 
that this prevents the peak at x < 1 from dominating for 
Pi(x), and as a result we obtain 4d space-time. Since the 
above argument is based only on the scaling behaviors 
and the branched polymer description, it is expected to 
be valid also in the D = 10 case. (While this paper was 
being revised, an analytic evidence for the dominance of 
4d space-time was also reported Q.) 

Our new approach to complex-action systems is based 
on the factorization property l|ll|) of distribution func- 
tions, which is quite general. As we discussed in Sec- 
tion H^21 h resolves the overlap problem completely. In 
a separate paper we will report on a test of the new 
method in a Random Matrix Theory for finite density 
QCD, where exact results in the thermodynamic limit 
are successfully obtained ^jjj. We hope that the 'factor- 
ization method' allows us to study interesting complex- 
action systems in various branches of physics. 
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FIG. 4: The function 4j>/4 (as) is plotted for N = 16,32,64, 128. A clear scaling behavior is observed. 
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